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The density matrix renormalization group (DMRG) algorithm was originally designed to efficiently 
compute the zero temperature or ground-state properties of one dimensional strongly correlated 
quantum systems. The development of the algorithm at finite temperature has been a topic of much 
interest, because of the usefulness of thermodynamics quantities in understanding the physics of 
condensed matter systems, and because of the increased complexity associated with efficiently com¬ 
puting temperature-dependent properties. The ancilla method is a DMRG technique that enables 
the computation of these thermodynamic quantities. In this paper, we review the ancilla method, 
and improve its performance by working on reduced Hilbert spaces and using canonical approaches. 
We furthermore explore its applicability beyond spins systems to t-J and Hubbard models. 
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I. INTRODUCTION 

The behavior of low dimensional strongly correlated 
electron systems at finite temperature has became a 
topic of great interest amongst the condensed matter 
physics community. Tuning the temperature can lead 
to interesting phenomena which cannot be fully under¬ 
stood in terms of ground state physics, such as the quan¬ 
tum critical regime,^ the transition to a spin incoher¬ 
ent Luttinger liquid,^ or the sudden emergence of a sin¬ 
gle spinon dispersion in XXZ-like spin-chain materials. 

On the experimental front, recent advances®^® have en¬ 
abled very precise measurements of both static and dy¬ 
namical correlation functions. These experimental re¬ 
sults are challenging the available theoretical and com¬ 
putational techniques,increasing demand to develop 
efficient and accurate numerical methods for the inves¬ 
tigation and eventual prediction of the thermodynamic 
quantities. 

In one dimension, the density matrix renormalization 
group^®^^® (DMRG) is the most powerful method to 
calculate ground state properties of strongly correlated 
systems. The DMRG has been successfully extended 
to treat real time dynamics^®^^® and systems at finite 
temperature. These extensions include transfer-matrix 
DMRG,^® minimally entangled typical thermal states or 
METTS,^'^^®® and the purification scheme^®’®®’^^^®® to be 
used in this paper. 

METTS was recently introduced by White in ref. 21, 
where he shows that these so-called minimally entangled 
typical thermal states efficiently represent the thermal 
properties of the system of interest. METTS has been 
applied to fermionic systems by one of us.®® Ref. 31 com¬ 
pares the purification scheme to the METTS method, 
concluding that METTS is, in general, computation¬ 
ally less efficient than purification. In particular, it is 
argued,®® that the additional statistical error source in¬ 
troduced by METTS sampling increases computational 
costs, especially at high temperatures. More recently. 


ref. 32 develops symmetry conserving METTS, improv¬ 
ing the efficiency of the method. Perhaps METTS and 
ancilla should be regarded as complementary. Given 
these recent developments in DMRG techniques for tem¬ 
perature dependence, the present paper revisits the an¬ 
cilla method, with an emphasis on improving its perfor¬ 
mance, and exploring its applicability beyond spins sys¬ 
tems to t-J and Hubbard models. 

The ancilla method was originally introduced within 
a thermoheld formalism by Umezawa®®, and later in 
DMRG®®’®'. It is designed to calculate the finite temper¬ 
ature quantum average of an observable, by definition a 
trace over all the eigenstates of the physical chain, and 
not accessible with standard DMRG. To this aim, the 
above quantity is expressed as a quantum average over a 
pure state— a thermal vacuum—in an enlarged Hilbert 
space. The Hilbert space is enlarged by adding an ancilla 
site to each site of the physical lattice, thus obtaining a 
two-leg ladder geometry. The ancilla’s degrees of free¬ 
dom are added in such a way that the reduced density 
matrix of the physical sites reproduces the thermal den¬ 
sity matrix. 

The purification scheme starts with an infinite temper¬ 
ature state, given by a product of maximally entangled 
states on each rung of the ladder. For spin chains,®® 
each ancilla can be chosen to have “opposite” quan¬ 
tum numbers to those of the corresponding physical site. 
This choice of ancillas is equivalent to applying a time- 
reversal transformation in the case of spin chains, or, 
more generally, a particle-hole transformation on the an¬ 
cillas, and has its profound motivation in giving the cor¬ 
rect Kubo-Martin-Schwinger relations®®^®® at finite tem¬ 
perature. The previously described choice has the addi¬ 
tional advantage of reducing the dimension of the Hilbert 
space of the total system composed of physical and an¬ 
cilla chains. For example, in the case of a Heisenberg 
model in the absence of a magnetic field, one is able to 
work in a subspace of the enlarged Hilbert space of physi¬ 
cal plus ancilla chains with z component of the total spin 
equal to zero. After the preparation of the initial state. 
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a finite temperature state is then obtained by evolving 
in imaginary time with the Hamiltonian of the physical 
chain. 

In this paper we show that the above approach, which 
we define as grand canonical for the physical chain, is not 
the most efficient for the DMRG numerical simulations: 
the preparation of the initial infinite temperature state 
and the subsequent imaginary time evolution can be fur¬ 
thermore restricted, without any additional truncation of 
the Hilbert space, to a subspace where the z component 
of the total spin of the physical chain is also conserved. 
This subspace has clearly smaller dimension than the one 
previously mentioned, the one that conserves only the 2 
component of the total (physical plus ancilla) system. 


II. THE ANCILLA METHOD FOR FINITE 
TEMPERATURE DMRG 

For convenience, we here summarize the ancilla ap¬ 
proach for finite temperature DMRG; refs. 16,37 have 
more detail. At /I = 0, that is, at infinite temperature, 
a maximally entangled state {ipi/d = 0)) between the 
physical sites of the chain and their ancillas is initially 
produced. A pure state in the enlarged system at finite 
temperature is then calculated by evolving |'!/:(/3 = 0)) 
in imaginary time with the Hamiltonian of the physi¬ 
cal sites, lipiP)) = = 0)). Given a generic 

observable O of the physical chain, the thermodynamic 
average can be calculated^®’^^^^® in the enlarged space 
using the standard zero temperature expression 


How can we impose the restriction outlined above? Or, 
more generally, how can we work in the canonical ensem¬ 
ble for the physical chain? To answer this question we 
will show how to engineer an “entangler” Hamiltonian 
that conserves a given set of quantum numbers of the 
physical chain, and generates a maximally entangled in¬ 
finite temperature state that is different from the state 
used in the grand canonical approach. Introduced by 
Feiguin and Fiete in ref. 25 for a t-J chain, it turns out 
that the canonical “entangler” Hamiltonian is non-local 
for the total system, but involves interactions between 
the rungs of the total system at all possible distances. 
In the thermodynamic limit, the canonical scheme gives 
the same results as the grand canonical one. Yet the 
performance of the DMRG simulation is remarkably im¬ 
proved, by up to one order of magnitude with respect to 
the grand canonical approach. 

Section HA introduces the grand canonical and the 
canonical purification schemes for the Heisenberg model. 
Section H B then extends the treatment to the case of the 
t-J model, and section HG to the Hubbard model. The 
Heisenberg model has only spin degrees of freedom, but 
the t-J and Hubbard models have also charge, leading to 
different “entangler” Hamiltonians. For each model, we 
first consider a local entangler, where neither the total 
number of electrons, nor the total number of spins of the 
physical chain are conserved; they are conserved only in 
average at finite temperature. This is the grand canonical 
ensemble. It needs a (temperature-dependent) chemical 
potential in order to keep constant the average number of 
electrons in the physical chain. A magnetic field is like¬ 
wise needed to keep constant the average z component 
of the total spin (as in the case of a spin chain). We then 
consider a non local entangler where the number elec¬ 
trons is conserved but not the z-component of the spin. 
Even in this case, we show that the z-component of the 
spin for the electrons in the physical chain is conserved 
in average during the temperature evolution; its value is 
zero in the absence of a magnetic field. Finally, we write 
down an entangler such that both the number of electrons 
and the z-component of the spin of the physical chain are 
conserved. 
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where the norm Z{f3) = (■!/)(/?) |i/)(/3)) represents the par¬ 
tition function at temperature T = 1/13. 

We present two different purification schemes for the 
initial infinite temperature state. The first scheme is re¬ 
ferred to as “grand canonical,” where one conserves only 
quantum numbers for the enlarged system given by the 
physical and the ancilla sites. In the second scheme, one 
conserves quantum numbers not just globally, that is, 
not just in the enlarged system, but in the physical and 
ancilla chains separately as well. 

In the next sections we will often discuss intensive ener¬ 
gies of the models considered. The definition of intensive 
energy is 


L mmw)) 
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where H is the Hamiltonian of the model under consid¬ 
eration and L is the system size. 


A. Heisenberg model 

1. Grand-canonical purification scheme 

Let us start considering the case of a spin chain de¬ 
scribed by the Heisenberg model 

L-2 

= (3) 

i=0 

with S = {S^, and where the Hilbert space of a 

single site is two-dimensional, having only two possible 
states, It) and |t). Recall briefly the purification scheme 
adopted in ref. 16, which in our classification turns out 
to have a grand-canonical character. For simplicity, first 
consider the case of two spins, accompanied by their an¬ 
cilla sites. The infinite temperature state is 

IV '2 spi„s(^ = 0)) = i(| n) + I tt)) G (I tt) + I tt)), (4) 
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where 0 is the direct product of states on the two com¬ 
posite sites, each given by a pair consisting of a physical 
site and its ancilla. In this paper, the first entry of the ket 
vector refers to the state of the physical site and the sec¬ 
ond entry to the ancilla site. For example, ket |ti), has 
a spin up in the physical site and a spin down in the cor¬ 
responding ancilla site. In the state Eq. (4), each ancilla 
has “opposite” quantum numbers with respect to those 
of the physical site. As long as the mapping between 
the states of a physical site and its ancilla is one-to-one, 
then any choice for the mapping of quantum numbers 
of the ancillas yields the correct thermodynamics for the 
physical chain. 

The choice outlined above (which is equivalent to ap¬ 
plying a time-reversal transformation^®) has the advan¬ 
tage of reducing the dimension of the Hilbert space of the 
total enlarged system. By construction, the state Eq. (4) 
lives in the subspace with spin 

of the total physical-plus-ancilla Hilbert space, where we 
have indicated with a(i) the ancilla site corresponding to 
the physical site i. Eor L sites, Eq. (4) generalizes to a 
so-called product of local maximally-entangled states 

|i/,(/3 = 0)) = ^ ^ Ictct), (5) 

i=0 cr=t4 


By diagonalizing a 4 x 4 matrix, where 4 is the dimension 
of the composite physical plus ancilla site, the entangler 
for the Heisenberg model in the grand-canonical scheme 
can be written as 

= ( 6 ) 

i=0 

We address the results coming from this entangler in 
the next section, where we compare the grand canonical 
and canonical purification schemes. 

2. Canonical purification scheme 

Let us consider the case of the Heisenberg model in 
the absence of a magnetic field. In this case, one would 
build a maximally entangled state with the property that 
SfoM) = 0. -S'ph.k) = 0, Sl^ \4>) = 0, such that one con¬ 
serves the spin of the physical and ancilla chains sepa¬ 
rately. 

For the simple case of two spins, a maximally entangled 
state with the above characteristics is 

\ 4 > 2 spi„s(/3 = 0)) = ^(in)01;t) + ut)®In))- (7) 


where Af is a normalization constant, f =4- ^-nd f =t. 
The above state is an eigenstate of the z-component of 
the total spin but it is not an eigenstate of the total 
spin of either the physical ancilla sites 

QZ — QZ 

^an. — ^a{i)- 

The above observation has important consequences for 
the efficiency of the finite temperature evolution. Indeed, 
within the purification scheme outlined above, one can 
study the thermodynamics of a generic spin chain Hamil¬ 
tonian in the presence of an external magnetic field. The 
average magnetization of the system can be tuned by 
changing the strength of the field at finite temperature. 
This possibility stems from the symmetry property of the 
initial infinite temperature state, which conserves neither 
^p\. nor Sl^,. 

A state like Eq. (5) is exponentially large, and would 
have to be truncated before time evolving it. Otherwise, 
one would have to work with a vector of size 2^. Eor the 
grand canonical scheme, one could build it by growing 
it slowly inside a traditional DMRG, and truncating it 
along the way. But due to the DMRG transformations, 
the book keeping of such a state would be unfeasible to 
perform. It is then clear that with traditional DMRG 
the entangler is needed for an efficient representation of 
the state. With an MPS approach, the infinite temper¬ 
ature state in Eq. (5) has bond dimension 2, and is not 
entangled in a global sense, because it is just a product 
state of singlets on each rung. How can we then gener¬ 
ate state Eq. (5) as a starting point for the temperature 
evolution? As proposed by Feiguin and Fiete in ref. 25, 
we find useful the notion of “entangler” Hamiltonian: a 
Hamiltonian having the state Eq. (5) as its ground state. 


Eor L sites 


IV'(/3 = 0))c 
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where M' is a normalization constant, =o) is the 

projector operator such that the z-component of the total 
spin of the physical (ancilla) chain is conserved and equal 
to zero: = 0))c = = 0))c = 0. 

For an L site chain, the subspace containing the max¬ 
imally entangled state in the canonical approach has di¬ 
mension (Cf gives the number of possible combina¬ 
tions of I objects on k places), a number clearly smaller 
than the analogous subspace dimension of the grand 
canonical case. As will be seen below, imposing the 
conservation of ^'pjj during the finite temperature evolu¬ 
tion in the enlarged Hilbert space substantially reduces 
the computational effort of the time dependent DMRG 
simulations.®® Conserving this symmetry can thus help 
improve DMRG implementations, because they usually 
consider global symmetries only. 

How can we, in practice, generate the maximally en¬ 
tangled “canonical” state Eq. (8)? By finding the correct 
“entangler” Hamiltonian, the one that gives Eq. (8) as its 
ground state. For the Heisenberg model, this “entangler” 
Hamiltonian is 


= (9) 

where fJ = S~^.y Notice that the above Hamiltonian 

has long-range interactions: each rung interacts with all 
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FIG. 1: (Color online) Intensive energies for a Heisenberg 
chain with J = — 1 and length L = 6 sites for the canonical 
and grand canonical approach, and comparing DMRG and 
ED, as indicated. 


the others. Section III discusses the entanglement impli¬ 
cations caused by long-range interactions in the DMRG 
calculations. A proof that the state (8) is the ground 
state of the Hamiltonian (9) is provided in appendix B. 

In matrix product operator notation,^® = 




/ 0 r| 
0/0 
or,/ 


for 0 < z < L — 1, = [/,0,rQ], and 


^R-i] 


0 

I 
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Fig. 1 shows the intensive energy defined as as a func¬ 
tion of the temperature for a ferromagnetic Heisenberg 
chain with L = 6 sites, and where we assume |J| = 1 
as unit of energy. Fig. 1 compares both purification 
schemes. It also compares DMRG with exact diago- 
nalization. The figure shows a perfect numerical agree¬ 
ment between DMRG and ED, but substantially differ¬ 
ent results between the grand canonical and canonical 
approaches; the intensive energies calculated within the 
canonical approach are systematically bigger than those 
calculated in the grand canonical scheme. Fig. 1 also 
shows that the results obtained with different approaches 
coincide for /? —>■ oo, that is, at zero temperature. In ap¬ 
pendix A, we prove that for a generic interacting quan¬ 
tum system the canonical (when chosen in a symmetry 
sector containing the ground state) and grand canonical 
average energies differ in general, but coincide at zero 


FIG. 2: (Golor online) Main panel: average intensive energies 
for the Heisenberg model (J = —1) at temperature T = 0.1 
as a function of the system size comparing the canonical and 
the grand canonical approaches (we consider up to L = 64 
sites). The Bethe ansatz yields an energy equal to —0.240219 
in the thermodynamic limit. Inset: difference between the 
intensive energies obtained in the two approaches, (AE) = 
(Ec) — (Eqc), as function of the temperature, for different 
system sizes: triangles (brown) indicate L = 6, diamonds 
(green) indicate L = 12, and crosses (cyan) indicate L = 24. 


temperature, and in the thermodynamic limit (where 
they coincide at any temperature). 

The DMRG implementation used throughout this pa¬ 
per is discussed in the supplemental material which can 
be found at https://web.ornl.gov/~gzl/papers/55/. 
The algorithm implemented for the time-dependent part 
of the DMRG is based on the Krylov space decompo¬ 
sition method, as explained in ref. 39 and references 
therein. An implementation can be found in ref. 40. In 
the time-evolution calculations, 600-800 states per block 
are kept. The results are well converged with this number 
of states in the entire interval of temperature investigated 
iPmin = 0, l3max/‘2. = 40, with A/3/2 = 0.1). 

Fig. 2 substantiates another analogy of the results ob¬ 
tained in the two purification schemes with those ob¬ 
tained in different ensembles. The main plot of fig. 2 
shows the average energies obtained in the canonical and 
grand canonical approaches as a function of the system 
size at low temperature, T = 0.1. The results are differ¬ 
ent for small system size, but they converge to the exact 
Bethe ansatz solution when extrapolated to the thermo¬ 
dynamic limit. The inset shows the difference between 
the average energies in the grand canonical and canonical 
approaches, which decreases at every temperature with 
increasing system size. On the other end, all the curves 
converge to the same value for /3 —>■ oo, as already ob¬ 
served in fig. 1. In the supplemental material, which can 
be found at https://web.ornl.gov/~gzl/papers/55/, 
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we have proved that this is valid for a generic quantum 
system when the thermodynamic limit is taken. 


B. t-J chains 

1. Grand canonical purification scheme 

In this section, we apply the theory discussed previ¬ 
ously to the case of a < — J model described by the stan¬ 
dard Hamiltonian 


L-2 

HtJ = -t (cl.crCi+l.cr + h.c) 

2—0,(7 

L-2 

+ j ^ - fV,iV,+i/4), (12) 

where Ni = X)<T=t ^ o-Ci.ir- The model is character¬ 
ized by a three-dimensional single-site Hilbert space with 
states empty |0), single occupied with spin up |t), and 
spin down j),). We start with the grand-canonical purifi¬ 
cation scheme. Following the same procedure introduced 
in the spin chains’ case, one can straightforwardly write 
down a product of local (along the rungs of the ladder) 
maximally entangled states: 


1^ = 0)) 


L-1 


2 = 0 


|0,0)-b 5] 


era- 




(13) 


We have verified (by diagonalizing a 9 x 9 matrix, where 
9 is the dimension of a composite physical-plus-ancilla 
site) that this state can be generated by calculating the 
ground state of the entangler Hamiltonian 

L-l 

H^oc = E (14) 

2=0 

where A| = (4,wl(dd “ observed in 

the case of the spin chains, even though Eq. (13) con¬ 
serves the total spin in the enlarged physical plus ancilla 
combined system, it does not conserve these quantities 
separately for the physical chain or for the ancilla chain. 
In order to get the thermodynamic properties for the 
physical chain at finite temperature in the grand canon¬ 
ical ensemble, one must add a chemical potential term /r 
during the imaginary time evolution. To keep an average 
of N electrons in the physical chain, one needs to solve 
the equation 

(V'(/3,Ai)|iVph.|V’(/3,M)) = A', (15) 

for each temperature T = 1//3, where |'!/'(/3, m)) = 
g-/3/2[//tj-/jArph.] with HtJ being the standard 

t — J Hamiltonian and iVph. = Yli -^z,ph- For this reason, 
when studying the thermodynamic properties of the t-J 
model at fixed density and zero magnetic field, the grand 


canonical scheme outlined above has a clear disadvantage 
in terms of computational cost. Therefore, results in the 
grand canonical scheme are calculated only with exact 
diagonalization for a small system size. We have verified 
that the average spin of the physical chain is zero at any 
temperature, ('0(/3,/r)|5'pij |'(/'(/3, A^)) = 0, with being a 
solution of Eq. (15). 


2. Canonical purification scheme 


In this section, we address the “canonical” purifica¬ 
tion scheme for the t — J model. We present two dif¬ 
ferent schemes which take into account the charge and 
the spin symmetries of the physical chain. We first re¬ 
view a canonical purification method where one only em¬ 
ploys the charge conservation in the physical chain. This 
scheme has already been treated by Feiguin and Fiete 
in ref. 25 in the context of a spin-incoherent Luttinger 
liquid. We consider a t-J chain with L sites, and a even 
number of electrons N, such that = N/2. No¬ 

tice that in the total system given by physical and ancilla 
sites one has 2N electrons. Up to a constant, we here 
provide the expression 


IV’(/3 — 0))ci — P(Npi,-N) 


L-l 


2—0 


(16) 

for the “canonical” maximally entangled state in the L 
sites case, where L^ph =n} is a projector operator (dif¬ 
ferent from the projector that appears in Eq. (8)) such 
that the total number of electrons of the physical (an¬ 
cilla) chain is conserved, Aph,|(/^(/? = 0))ci = ./Van.|'0(/? = 
0))ci = N. We emphasize that the above state does not 
conserve the the z-component of the total spin of the 
physical (ancilla) chain, but it has the property that, at 
any temperature, ci(V'(/5)|5'ph.lV’(^))ci = 0. 


We now prove that the “canonical” entangler Hamil¬ 
tonian 


iJ^7-^ = -EAlA,+h.c. (17) 

proposed in ref. 25 has the property that its ground-state 
is Eq. (16). We begin by observing that the Hamiltonian 
(17) conserves the number of electrons N in the phys¬ 
ical and ancilla chains separately, and assume that the 
combined physical plus ancilla system has = 0, so 
that the number of up and down electrons are equal in 
the combined system. Let us divide the full Hilbert space 
basis into “good states” where all physical sites and their 
ancillas are correctly paired, and “bad states” where at 
least one physical site is not correctly paired to its ancilla; 
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these disjoint sets are then 

Sc ={|(^);basis \(p) with all physical and ancilla 
sites correctly paired} 

Sb ={|(^); basis |0) with at least one physical site 
with ancilla incorrectly paired}. 


Note that if \(j)) S Sc then |0) is a term in Eq. (16). We 
shall prove that (i) {(j)'\H^'^{(j)) = 0 if \(j)') G Sc and 
|(()) S Sb, so that the Hamiltonian matrix blocks 

into at least two blocks: states in Sc and states in iS^. 
We shall furthermore prove (ii) that the ground state of 
is in the block Sc as opposed to the block Sb, and 
(iii) that the ground state of the block 16) for 

all !<).), If) is Eq. (16). 

To prove (i) it suffices to think of Eq. (17) as a tight- 
binding Hamiltonian of “singlets” of physical sites and 
ancillas built along the rungs, singlets that are hopping 
along the ladder structure. Eq. (17) cannot connect 
states in Sc with those in iS^. 

To prove (ii) we first note that if \(j)) G Sb, then at 
least one site has occupation |cr, 0), |0, cr) or |cr, cr). Elec¬ 
trons in sites like these cannot “move” by the action of 
Eq. (17). The subspace Sb then blocks into even smaller 
subspaces, all of size smaller than the set Sc- Within 
these subspaces, we have hopping-like matrices yielding 
lowest energies larger than that of the block Sc- 

To prove (iii), let us call H the block of matrix Eq. (17) 
in the Sc subspace. The Hamiltonian Eq. (17) conserves 
the number of electrons in the physical system, but not 
its total spin, so that it can be written in a tridiagonal 
block form. The diagonal blocks have rank 
each, and are characterized by the number of up elec¬ 
trons. The blocks in the first diagonal above (below) 
the main diagonal connect states with and -G 1 
{N^ — 1) electrons and are given by m x n rectangular 
matrices with m = and n = 

(r, — 

Now, the matrix elements in the blocks described 
above are either 0 or -1 for the diagonal blocks while for 
the off-diagonal ones are 0 and -1-1. Each row has exactly 
2N{L — N) non-zero entries. So does each column. Then 
its lowest eigenvalue (by taking into account the factor 
1/2 in the definition of the A operators) is —N{L — N), 
with eigenstate (1,1, • • • , 1), that is, Eq. (16). 

Because the t — J chain has charge and spin degrees of 
freedom, there is an alternative and more efficient way of 
performing the canonical purification scheme: Consider 
the maximally entangled infinite temperature state 

' L-l 

l^(/3 = 0 ))c 2 = C’/-«ph,=«') (^(|0,0)-f- ^ \ad)) , 

^ L i=0 <T=td 

(19) 

where the projector operator is such that the total num¬ 
ber of electrons and the total z-component of the physical 



T 

FIG. 3: (Color online) Intensive energies for a t-J chain of L = 
6 sites, N — 4: electrons, with J = 0.1, for the canonical and 
grand canonical approaches (Cl and C2, see main text), and 
comparing DMRG against ED, as indicated. Inset: Chemical 
potential as a function of the temperature, calculated in the 
grand canonical approach so that the average density is always 
N = 4 electrons. 


(ancilla) chain are conserved, 

TVph. |^(/3 = 0))c2 =N\i,il 3 = 0))c2 (20) 

N^n. mP = 0))C2 =NmP = 0))C2 (21) 

m = 0))c2 mp = 0))c2 = 0. (22) 

With a procedure similar to the one just outlined, it is 
possible to show that an “entangler” Hamiltonian yield¬ 
ing Eq. (19) as ground state is 

~ ^ +^- C-- (23) 

The ground state energy of the above Hamiltonian is 
equal to —N{L — N), coinciding with the result for 
the Cl Hamiltonian. Notice that, for an L site chain, 
the subspace containing the maximally entangled state 
in the canonical approach Cl has dimension Dqi = 
'which is evidently larger than the 

analogous subspace dimension Dc 2 = 
the case C2, where one has also = N/2. If 

one imposes the conservation of Aph. and during 
the finite temperature evolution in the enlarged Hilbert 
space, this remarkably reduces the computational effort 
and increases the efficiency of the time dependent DMRG 
simulations. In our typical runs, when this purification 
scheme is adopted, the computational time needed for 
obtaining the thermodynamic properties is reduced by a 
factor of two with respect to the canonical purification 
Cl. 
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Fig. 3 shows the average intensive energy as a func¬ 
tion of the temperature for a t-J chain with L = 6 sites, 

= 4 electrons, and J = 0.1, (t = 1 is assumed as unit 
of energy) comparing canonical and grand canonical ap¬ 
proaches. One observes a perfect numerical agreement 
between DMRG and ED for both canonical approaches 
Cl and C2. For the small system size considered, the 
results for the average energies differ (by less than 1%) 
at temperature larger than zero, while converging to the 
ground state value in the limit of large /3. We have ver¬ 
ified that, even for a system size of L = 12 sites, the 
numerical results obtained in the Cl and C2 approaches 
are coincident within the numerical precision of our runs. 

The results obtained in the grand canonical approach 
were calculated with ED by performing a procedure sim¬ 
ilar to a standard Maxwell construction. We first cal¬ 
culated the total electronic density as a function of the 
temperature, for different chemical potential values /r. 
By solving an equation similar to Eq. (15), (Ntot) = 4, 
we have then extracted the chemical potential curve at 
constant density as function of the temperature, which 
is reported in the inset of fig. 3. Finally, in the main 
plot of fig. 3, the blue crosses represent the intensive en¬ 
ergy as a function of the temperature. Even for small 
system sizes, DMRC calculations in the grand canonical 
approach are already computationally expensive. We can 
then conclude that the canonical purification C2 should 
be the computationally preferred choice for the study of 
the thermodynamic properties of the t — J model. (Yet 
the disadvantages of the canonical are discussed in sec¬ 
tion III.) 

C. Hubbard chains 


Eq. (13) is that the double occupied state on a physical 
site is “mapped” to the same state on its ancilla. By 
diagonalizing a 16 x 16 matrix, where 16 is the dimension 
of a composite physical-plus-ancilla site, the state above 
can be generated by calculating the ground state of the 
entangler Hamiltonian 

L-l 

^Hubbard ^ ^ +h.C.), (26) 

i = 0.dr=t4 

where Pf = |1 - N,^g. - fVa(z),cri¬ 
state Eq. (25) conserves the total spin in the enlarged 
chain, but does not do so for the physical chain. A chem¬ 
ical potential /i must therefore be added in the Hamilto¬ 
nian used during the temperature evolution, and mutatis 
mutandis, the same procedure outlined in the previous 
section for the t — J case must be followed. The average 
total spin of the physical chain is zero at any temper¬ 
ature, that is, fi)) = 0, with /i being 

a solution of an equation similar to Eq. (15). Because 
of the computational cost of the procedure, the results 
in the grand canonical scheme are presented using only 
exact diagonalization on a small system size. 

2. Canonical purification scheme 

We now address the “canonical” purification schemes 
for the Hubbard model. We consider a chain with L sites, 
filling N/L, and a total even number of electrons N, such 
that N^ = N^= N/2. 

Up to a constant, the “canonical” maximally entangled 
state of type Cl at infinite temperature is 


1. Grand canonical purification scheme 

In this section, we apply the purification schemes in¬ 
troduced in the previous section to a one dimensional 
Hubbard model described by the standard Hamiltonian 

L-2 L-l 

Hnub = -^ {c\^„Ci+i,c+h.c)+U'^Ni^-^Ni^^. (24) 

Z—0,(7 z—0 

The model is characterized by a local four-dimensional 
Hilbert space with states empty |0), single occupied |t), 
II) and double occupied |tl)- 
We start by considering the grand-canonical purifica¬ 
tion scheme. Following the same procedure used for the 
t — J case, one can straightforwardly write (up to a con¬ 
stant) 


|^)ci = ^(iVph-iV) 


"L-l 

(|0,0) +1 ti, ti) + ^ 1^^)) 


. z=0 


o-=t,l 


(27) 

where P(iVph=Af) is a projector operator such that the 
total number of electrons of the physical (ancilla) chain 
is conserved. 


-^ph. |f/’,S=o)ci = A^an.|V'/3=o)ci = -Y|V'/3=o)ci• (28) 

Even though this state conserves charge in the physical 
chain, it does not conserve the z-component of the spin 
of the physical chain. Yet, ci(l/’(/3)|*S'ph.|t/'(/3))ci = 0. 

As shown below, the canonical purification just dis¬ 
cussed is not the most efficient for a DMRC implementa¬ 
tion. For this reason, in the rest of the section we focus 
on the canonical scheme of type C2: the maximally en¬ 
tangled infinite temperature state 


L-l 


i=Q 


io, 0) +1 tiiti) + 


(25) 




for the infinite temperature state in the grand canonical 
approach. The only difference with respect to the state 


|V’)C2 = -Pl«ph.=«'| 


L-l 


(|o, 0) -I- 111,14-) + 1®’'^)) 


i=0 


a=t,l 


(29) 

where the projector operator is such that the total num¬ 
ber of electrons and the total z-component of the physical 
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FIG. 4: (Color online) Intensive energies for the Hubbard 
model, U = 10, and density 2/3, as a function of the system 
size. Squares indicate data obtained with DMRG using a 
maximum of m = 1000 and a truncation error of 10“®. The 
dashed (magenta) lines indicate the energy calculated with 
thermodynamic Bethe ansatz. Left panel shows T = 1.25, 
right panel T = 0.25. 


(ancilla) chain are conserved, 

^^ph.|lA/3=o)c2 =./Van.|l/’/3=o)c2 = N\tpp=o)c 2 , , . 

^p\.|lA/3=o)c2=5f„.l#=o)c2 = 0. ^ ^ 

With a procedure similar to that outlined in section 
IIB 2, it is possible to show (see appendix C) that an 
“entangler” Hamiltonian giving Eq. (29) as ground state 
is 

= - Y. ( 31 ) 

where o- = Ci.sCi^rrPi ^ where is given below 
Eq. (26). ’ 

For an L site chain, the subspace containing the max¬ 
imally entangled state in the canonical approach Cl has 
dimension Dqi = which is larger than 

the analogous subspace dimension Dc 2 = Cn/ 2 ^n /2 
the case C2; both cases have = N/2. (Notice 

the difference with the corresponding expressions pro¬ 
vided for the t — J model in the previous section.) If 
one imposes the conservation of A^ph. and S'pj, during 
the finite temperature evolution in the enlarged Hilbert 
space, one can thus reduce (in our typical runs by a fac¬ 
tor of two or more) the computational time needed for 
obtaining the thermodynamic properties. We therefore 
recommend that the canonical approach C2 be the pre¬ 
ferred purification scheme when finite temperature static 
and even dynamic properties are calculated with DMRG. 

Figure 4 shows the intensive energies as a function of 
the system size for the canonical approach C2 for U = 
10 (t = 1 is assumed as unit of energy) and filling 2/3. 
The left panel shows data at T = 1.25, while the right 
panel at T = 0.25. Using m = 1000 and imposing a 
truncation error not bigger than 10“® in the DMRG runs, 
the results get closer to the thermodynamic Bethe ansatz 



FIG. 5: (Golor online) Intensive energies for the Hubbard 
model, U = 10, and density 2/3, as a function of the temper¬ 
ature for a L = 6 chain. Shown are results within the grand 
canonical, calculated using ED, blue crosses; canonical Cl, 
calculated using ED, black squares; canonical C2 calculated 
using ED, green triangles; and canonical C2 obtained with 
DMRG, dark green empty circles. The ground state solution 
is also indicated with a magenta dashed line. 


results as the system size increases. In this paper, we 
have considered chains of length up to L = 60. 

Fig. 5 shows a comparison between the intensive en¬ 
ergies obtained in the canonical and grand canonical ap¬ 
proaches for a Hubbard chain with L = 6 sites, = 4 
electrons, and 17 = 10 as a function of temperature. Be¬ 
cause the system size is small, the numerical results differ 
in the three approaches, even though they converge to the 
exact ground state solution in the limit of /3 —> oo. Re¬ 
markably, even for L = 6, the results of the G1 and C2 
approaches differ by less than 5% on the entire /3 interval 
studied. 

For the canonical purification scheme G2, Fig. 5 also 
compares DMRG results against those obtained with 
ED. The agreement is close to the numerical precision. 
Finally, the computationally expensive results with the 
grand canonical approach are also shown. These results 
are calculated by performing a procedure similar to a 
standard Maxwell construction. 


III. LONG RANGE INTERACTIONS IN 
CANONICAL ENTANGLERS 

In this section, we discuss the properties of the entan¬ 
gler Hamiltonians in the canonical purification scheme. 
As already mentioned in the previous sections, the en¬ 
tangler Hamiltonians have long range interactions, with 
connections between sites at all possible distances. The 
resulting entanglement growth makes it difficult to com¬ 
pute the ground state of such Hamiltonians with the 
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CPU time (seconds) 


FIG. 6: (Color online) CPU time of typical runs in the grand 
canonical and grand canonical approaches for the Heisenberg 
model at /? = 0.8, J = — 1, and different system sizes. The 
black bars indicate the CPU time needed to solve the entan- 
gler Hamiltonian (for the GC see Eq. (6), for C see Eq.( 9)), 
while the red bars indicate the same quantity for the temper¬ 
ature evolution. In the top panel {L = 24) and middle panel 
(L = 48) m = 400 DMRG states are used for the temperature 
evolution, while in the bottom panel {L — 64) m = 600. The 
CPU time for GC with L=64 has been truncated (the actual 
value is 595,000), because it is larger than the maximum CPU 
time shown in this scale. 

DMRG. Therefore, one would think that the grand 
canonical approach, where the local entangler Hamiltoni¬ 
ans are used, would be more efficient for the calculation 
of the thermodynamic quantities. But the canonical pu¬ 
rification scheme is computationally much more efficient. 

Fig. 6 explains why. It shows the CPU times of the 
two purification schemes for the Heisenberg model. CPU 
times are divided in two parts: the calculation of the 
ground state of the entangler Hamiltonian (black bars), 
and temperature evolution up to the desired inverse tem¬ 
perature /? = 0.8 (gray (red) bars), starting from the 
ground state of the entangler Hamiltonian. The compu¬ 
tational time needed to calculate the ground state of the 
entangler in the grand canonical approach is too small to 
be visible in the figure. The CPU time needed for the 
canonical entangler Hamiltonian is larger and visible as 
black bars. The apparent disadvantage of the canonical 
purification versus the grand canonical scheme seems to 
increase if large system sizes are considered. For L = 64 
sites and m = 200, the CPU time needed in the canonical 
scheme for the ground state calculation is about three or¬ 
ders of magnitude longer than that needed in the grand 
canonical approach. 

As indicated by the gray (red) bars in Fig. 6, the CPU 
time of the temperature evolution of the canonical pu¬ 
rification is remarkably smaller than that of the grand 
canonical, because the temperature evolution takes place 
in a Hilbert space of reduced dimensions. (This results 



FIG. 7: (Color online) Difference between the ground state 
energy of the entangler Hamiltonian Eq. (31) for the Hubbard 
model at density 2/3 and the exact ground state energy as 
a function of the number of DMRG sweeps. A chain with 
L = 60 sites is considered. Increasing the number of DMRG 
states kept from m = 200 to m = 400 allows the results to 
converge to the exact ground state energy. 


from the conservation of as stated in sec. H A.) Given 
that the total CPU time is equal to the sum of the CPU 
time needed for the entangler plus the CPU time needed 
for the temperature evolution, we can conclude that the 
canonical purification scheme is faster than the grand 
canonical approach—by one order of magnitude or more. 

Yet it is important to make sure that ground states of 
canonical entanglers are converged, given the aforemen¬ 
tioned difficulties due to the long range interactions. To 
that aim. Fig. 7 is an example showing that the canoni¬ 
cal entanglers used throughout this paper are converged. 
Fig. 7 shows the ground state energy of the canonical en¬ 
tangler C2 for the Hubbard model as a function of the 
number of DMRG sweeps, for a chain of L = 60 sites 
with U = 10 and density 2/3. Several sweeps and two 
values of m were used to converge to the exact ground 
state. 

We end this section by providing the ground state ener¬ 
gies of the entangler Hamiltonians. Let L be the number 
of physical sites of the system. As shown in Appendix B, 
the ground state energy of the Heisenberg entangler Eq. 9 
is -L^/4. Let N be the number of electrons in the phys¬ 
ical system. Section H.B.2 shows that the ground state 
energy of the t-J entangler Eq. (17) is —N{L — N), which 
is equal to that of Eq. (23). Appendix C proves that the 
ground state energy of the Hubbard entangler Eq. (31) 
is -N{L-N/2). 
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IV. SUMMARY AND CONCLUSIONS 

In summary, we have improved the efficiency of the 
ancilla method for finite temperature DMRG by em¬ 
ploying the inherent symmetries of the physical system 
in consideration. We have designed different entangler 
Hamiltonians to obtain infinite temperature states to use 
as starting states for the temperature evolution. Ta¬ 
ble I reports the purification schemes adopted for each 
model considered, and serves as an index to the equations 
obtained. The supplemental material can be found at 
https://web.ornl.gov/~gzl/papers/55/. It provides 
a pointer to the full open source code, input decks and 
additional computational details. 

As a main result, we have derived entangler Hamiltoni¬ 
ans for the canonical purification scheme of spin chains, 
the t — J model and the Hubbard model. The present 
work codifies an efficient ancilla method for spin chains 
and fermionic systems, because (i) canonical purification 
is the most efficient for obtaining the thermodynamic 
properties of the physical system in consideration, and 
(ii) entangler Hamiltonians are needed to obtain workable 
representations of infinite temperature states. Moreover, 
the efficiency brought about by the use of symmetries 
does not compromise accuracy: the ancilla method is as 
accurate as originally proposed. 

Due to their long range interactions, canonical entan- 
glers appear computationally costlier than grand canon¬ 
ical ones. But grand canonical entanglers have to be 
simulated on larger Hilbert spaces. Overall, canonical 
entanglers turn out to be computationally much more ef¬ 
ficient, as was shown in section HI. 

In the Hubbard model away from half-filling, we have 
verified that the efficiency gain of the canonical scheme 
overcomes the extra cost due to the need for more states 
in the temperature evolution for small and medium-size 
systems (up to L = 60). For very large system sizes 
and at half filling, the grand canonical approach remains 
an alternative method for the calculation of the thermo¬ 
dynamic properties'^. Yet the grand canonical scheme 
away from half-hlling requires computational runs to ad- 


Nt' 

^.ph. 

Heisenberg 

t-J Model 

Hubbard 

No 

No 


GC, (13), (14) 

GC, (25), (26) 

Yes 

No 

GC, (5), (6) 

Cl, (13), (17) 

Cl, (27) 

Yes 

Yes 

C, (8), (9) 

C2, (19), (23) 

C2, (29), (31) 


TABLE I: For each row, the conservation of the number of 
particles AJ’*' or number of particles AJ’*' and the z compo¬ 
nent of the spin of the physical system is indicated by 
“Yes” or “No.” For each model, the labels grand canonical 
(GC), canonical type 1 (Cl) or type 2 (C2) are defined. The 
numbers refer to equations in this paper: the first number 
is the equation defining the infinite temperature state, the 
second number (if present) is the equation defining the corre¬ 
sponding entangler Hamiltonian. For the Heisenberg model 
the A^'’ column should be ignored. 


just the chemical potential at every temperature! 

The purification schemes proposed are general and can 
be applied to both more complex one dimensional models 
and to geometries beyond chains. They apply to more 
general models because the Hamiltonian of the system 
comes into play only in the temperature evolution. They 
apply to more general geometries because the ancillas can 
be thought of as an extra orbital for each physical one. 

Concerning the Hubbard and t-J models, we believe 
that the canonical approach C2 should be the preferred 
purification scheme when finite temperature static and 
even dynamic properties are calculated with DMRG. For 
instance, in ref. 28, the entanglement growth character¬ 
izing the purification scheme has been reduced by time 
evolving the auxiliary degrees of freedom backward in 
time, when a combination of finite temperature and time 
dependent DMRG is needed for the calculation of spec¬ 
tral properties. We think that a natural step would be to 
combine the suggested improvement with the canonical 
purification proposed in this paper. 
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Appendix A: Canonical vs. grand canonical 
ensemble 

The purpose of this appendix is three-fold. First, to 
show that, for finite systems at finite temperature, the 
canonical and grand canonical ensembles do not, in gen¬ 
eral, yield the same averages. Second, to show that, at 
T = 0, the canonical and grand canonical ensembles 
agree for any Hamiltonian. And third, that they also 
agree in the thermodynamic limit for any Hamiltonian. 
Even though these results are known, we prove them be¬ 
low for completeness. 

Observation 1. At finite temperature, there exists at 
least one finite dimensional Hamiltonian with canonical 
average energy different than its grand canonical average 
energy. 

Proof. We construct such a Hamiltonian by considering 
2 levels, and a one-site Hilbert space including the states, 
empty, level 1 occupied only, level 2 occupied only, and 
both levels occupied. Let 

H = e{hih 2 )Vhih 2 , (Al) 
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where fii acting on the basis states multiplies it by the 
level I occupation of that state. Let e > 0 and V > 0. 
The average energy in the canonical ensemble with ex¬ 
actly = 1 particle in the system is {E)c = e. The aver¬ 
age energy in the grand-canonical ensemble with density 
(N) = 1 is given by {E)gc = e+ i+Jfil/2 , which is greater 
than e at finite temperature. □ 

Observation 2. For any finite dimensional Hamiltonian 
with convex energies, at zero temperature, the average 
energy in the canonical ensemble is the same as the one 
in the grand canonical ensemble. 

Proof. Let us consider a system with M sites and target 
density Nt. The number M includes sites, orbitals and 
spin, such that the maximum number of electrons that 
the system can hold is M. Let the eigenvalues of the 
Hamiltonian be E^,, where the index n' runs only over 
states of Fock sector J'jv of N particles. Let iJ^in = 
min„/gjrjy E^,. At T = 0, the average energy in the 
canonical ensemble with density Nx is E^f^. In the grand 
canonical (GC) ensemble we impose {N)gc,t=o = Nt, 
and thus the chemical potential p, is obtained by impos¬ 
ing that the N that minimizes F/v = — p.N be Nt. 

Because the energies are convex, such a /i is unique. For 
this fi, we obtain (F)gc.t=o = which coincides 

with the canonical result. □ 

Observation 3. For any Hamiltonian with convex en¬ 
ergies and an extensive canonical partition function, at 
finite temperature, and in the thermodynamic limit, the 
average energy in the canonical ensemble is the same as 
the one in the grand canonical ensemble. 

Proof. The proof is given in the supplemental material 
at https://web.ornl.gov/~gzl/papers/55/. □ 


Appendix B: Canonical entangler for the Heisenberg 

model 

We now prove that the “canonical” entangler Hamilto¬ 
nian for the Heisenberg model, Eq. (9), has the property 
that its ground-state is Eq. (8). We begin by observ¬ 
ing that the Hamiltonian (9) conserves the z-component 
of the total spin in the physical and ancilla chains sepa¬ 
rately, and assume that the combined physical plus an¬ 
cilla system has Sz = 0, so that the number of up and 
down spins are equal in the combined system. Let 

7g ={|(/)); basis |(/>) with all physical and ancilla 
sites correctly paired} 

Tb ={|(/'); basis |(/>) with at least one physical site 
with ancilla incorrectly paired}, 

(Bl) 

where Tg is the set of states in which each site is cor¬ 
rectly paired with ancilla according to the scheme 


}—)■}, and Tb the set where at least one site is incorrectly 
paired. 

It is easy to prove (i) that = 0 if |^') S 7 g 

and \(j)) € Tb, so that the Hamiltonian matrix 
blocks into at least two blocks: states in Tg and states 
in Tfi. To this aim, it is sufficient to think of spins as 
mapped to hard-core bosons, | t) = |1) I 4-) = |0)- One 
can imagine Eq. (9) as a tight-binding Hamiltonian of 
“singlets” of physical sites and ancillas built along the 
rungs, singlets that are exchanging positions of parti¬ 
cles and holes. Because Eq. (9) conserves the number 
of bosons in the physical and ancilla chain separately, it 
cannot, therefore, connect states in Tg (number of bosons 
and holes is equal to L/2 in the physical chain) with 
those in 7 b. In fact, in the subspace Tb at least one 
site of the chain has, in the hard-core bosonic represen¬ 
tation, occupation |0,0) (| i,}.)), |1,1) (| bosons in 

sites like these cannot “move” by the action of Eq. (9). 
As in the other cases, the subspace Tb then blocks into 
even smaller subspaces, all of size smaller than the set 
Tg- Within these subspaces we have hopping-like matri¬ 
ces yielding lowest energies which are larger than those 
in the set Tg- 

To prove (hi), let us call H the block of matrix Eq. (9) 
in the Tg subspace. The rank of H is (7^/2 
matrix elements are either 0 or -1. Each row has exactly 
L^/4 non-zero entries. So does each column. Then its 
lowest eigenvalue is —L^/4, with eigenstate (1,1, • • • , 1), 
that is, Eq. (8). 


Appendix C: Canonical entangler for the Hubbard 
model 

We now prove that the “canonical” entangler Hamilto¬ 
nian of type C2 for the Hubbard model, Eq. (31), has the 
property that its ground-state is Eq. (29). We begin by 
observing that the Hamiltonian (31) conserves the num¬ 
ber of electrons and z-component of the total spin in the 
physical and ancilla chains separately, and assume that 
the combined physical plus ancilla system has Sz = 0, so 
that the number of up and down electrons are equal in 
the combined system. As done in Eq. 18, let 

bVc ={|(/)); basis \(j)) with all physical and ancilla 
sites correctly paired} 

kVs ={|0);basis \(j)) with at least one physical site 
with ancilla incorrectly paired}, 

(Cl) 

where Wg is the set of states in which each site is cor¬ 
rectly paired with ancilla according to the scheme 0 —>■ 0, 
f—T'}, and Wb the set where at least one 

site is incorrectly paired. 

Even in this case, it is easy to prove (i) that 
{(j)'\H^2^\(j)) = 0 if \(j)') G Wg and \(j)) G Wb, so that 
the Hamiltonian matrix Hffff^ blocks into at least two 
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blocks: states in Wg and states in Wb- Indeed, it is 
sufficient to imagine Eq. (31) as a tight-binding Hamilto¬ 
nian of generalized “singlets’ of physical sites and ancillas 
built along the rungs 

I'i/’rung singlet) = |0, 0) -|- | + I t j 4-) + I J-it)- (C2) 

Because Eq. (31) conserves the number of electrons in 
the physical and ancilla chain separately, it cannot, there¬ 
fore, connect states in Wg with those in Ws. In fact, in 
the subspace Wb at least one site of the chain has occu¬ 
pation |(T,0), |0,cr), I tijO-), |cr,ti) or Electrons in 

sites like these cannot “move” by the action of Eq. (31). 


As in the other cases, the subspace Ws then blocks into 
even smaller subspaces, all of size smaller than the set 
Wg- Within these subspaces we have hopping-like ma¬ 
trices yielding lowest energies which are larger than those 
in the set Wg- 

To prove (iii), let us call H the block of matrix Eq. (31) 
in the Sq subspace. It can be written as a direct product 
H = ® H^. each has rank and 

matrix elements either 0 or -1. Each row has exactly 
N{L — NI2)I2 non-zero entries. So does each column. 
Then the lowest eigenvalue of H is —N{L — N/2), with 
eigenstate (1,1, • • • , 1), that is, Eq. (29). 
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